clear;
input_layer_size = 400; % 20x20 Input Images of Digits
hidden_layer_size = 25; % 25 hidden units
num_labels = 10; % 10 labels, from 1 to 10
load('ex4data1.mat');
m = size(X, 1);
sel = randperm(size(X, 1));
sel = sel(1:100);

displayData(X(sel, :));
load('ex4weights.mat');
nn_params = [Theta1(:); Theta2(:)];

lambda = 0;

J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, ...
    num_labels, X, y, lambda);

fprintf(['Cost at parameters (loaded from ex4weights):%f '...
    '\n(this value should be about 0.287629)\n'], J);

lambda = 1;

J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, ...
    num_labels, X, y, lambda);

fprintf(['Cost at parameters (loaded from ex4weights):%f '...
    '\n(this value should be about 0.383770)\n'], J);

initial_Theta1 = randInitializeWeights(input_layer_size, hidden_layer_size);
initial_Theta2 = randInitializeWeights(hidden_layer_size, num_labels);

initial_nn_params = [initial_Theta1(:); initial_Theta2(:)];

checkNNGradients;
lambda = 3;
checkNNGradients(lambda);

lambda = 3;
debug_J = nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda);
fprintf(['\n\nCost at (fixed) debugging parameters (w/ lambda = %f): %f ' ...
    '\n(for lambda = 3, this value should be about 0.576051)\n\n'], lambda, debug_J);

options = optimset('MaxIter', 50);
lambda = 1;
costFunction = @(p) nnCostFunction(p, input_layer_size, hidden_layer_size, num_labels, X, y, lambda);
[nn_params, cost] = fmincg(costFunction, initial_nn_params, options);
Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), hidden_layer_size, (input_layer_size + 1));
Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), num_labels, (hidden_layer_size + 1));

displayData(Theta1(:, 2:end));

pred = predict(Theta1, Theta2, X);
fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);